#pragma once

#include "sycl_devices.hpp"

/**
 * @brief  Initialize Fluid states espically primitive quantity;
 * @return void
 */
extern void InitialStatesKernel(int i, int j, int k, Block bl, IniShape ini, MaterialProperty material, Thermal thermal,
                                real_t *u, real_t *v, real_t *w, real_t *rho, real_t *p, real_t *_y, real_t *T)
{
    MARCO_DOMAIN_GHOST();
    if (i >= Xmax)
        return;
    if (j >= Ymax)
        return;
    if (k >= Zmax)
        return;

    int id = Xmax * Ymax * k + Xmax * j + i;
    real_t dx = bl.dx, dy = bl.dy, dz = bl.dz;
    real_t x = bl.DimX ? (i - Bwidth_X + bl.myMpiPos_x * (Xmax - Bwidth_X - Bwidth_X)) * dx + _DF(0.5) * dx + bl.Domain_xmin : _DF(0.0);
    real_t y = bl.DimY ? (j - Bwidth_Y + bl.myMpiPos_y * (Ymax - Bwidth_Y - Bwidth_Y)) * dy + _DF(0.5) * dy + bl.Domain_ymin : _DF(0.0);
    real_t z = bl.DimZ ? (k - Bwidth_Z + bl.myMpiPos_z * (Zmax - Bwidth_Z - Bwidth_Z)) * dz + _DF(0.5) * dz + bl.Domain_zmin : _DF(0.0);

    rho[id] = _DF(0.0);
    p[id] = _DF(0.0);
    T[id] = _DF(0.0);
    u[id] = _DF(0.0);
    v[id] = _DF(0.0);
    w[id] = _DF(0.0);

    // // 1D multicomponent insert shock tube
    if (x < ini.blast_center_x) // 0.05) // ini.blast_center_x) // (i > 3)
    {
        T[id] = ini.blast_T_in;
        p[id] = ini.blast_pressure_in;
        // T[id] = 400, p[id] = 8000;
        u[id] = ini.blast_u_in;
        v[id] = ini.blast_v_in;
        w[id] = ini.blast_w_in;
    }
    else
    {
        T[id] = ini.blast_T_out;
        p[id] = ini.blast_pressure_out;
        // T[id] = 1200, p[id] = 80000;
        u[id] = ini.blast_u_out;
        v[id] = ini.blast_v_out;
        w[id] = ini.blast_w_out;
    }

#ifdef COP
    real_t *xi = &(_y[NUM_SPECIES * id]);
    for (size_t nn = 0; nn < NUM_SPECIES; nn++)
        xi[nn] = _DF(.0);

    // for 2D/3D shock-bubble interactive
    real_t dy_ = _DF(0.0), tmp = _DF(0.0);

    if (bl.DimX)
        tmp = (x - ini.cop_center_x) * (x - ini.cop_center_x), dy_ += tmp * ini._xa2;

    if (bl.DimY)
        tmp = (y - ini.cop_center_y) * (y - ini.cop_center_y), dy_ += tmp * ini._yb2;

    if (bl.DimZ) // for 3D shock-bubble interactive
        tmp = (z - ini.cop_center_z) * (z - ini.cop_center_z), dy_ += tmp * ini._zc2;

    // // Ini bubble
    dy_ = sycl::sqrt(dy_) - _DF(1.0); // not actually the same as that in Ref: https://doi.org/10.1016/j.combustflame.2022.112085
    real_t xrest = _DF(1.0), ff = _DF(1.0e-4), dd = _DF(0.5) * (xrest - _DF(2.0) * ff);
    xi[NUM_SPECIES - 1] = _DF(0.0);
    xi[NUM_SPECIES - 2] = dd * (sycl::tanh(dy_ * ini.C)) + _DF(0.5); // increase ini.C for a sharper boundary //[-1,1]--0.5(1-ff)*[-1,1]+0.5-->[ff,1-ff]
#if defined(SBI_WITHOUT_FUEL)
    xi[0] = _DF(0.51) * (xrest - xi[NUM_SPECIES - 2]);               // O2
    xi[NUM_SPECIES - 3] = _DF(0.49) * (xrest - xi[NUM_SPECIES - 2]); // Xe
#else
    xi[0] = _DF(0.29) * (xrest - xi[NUM_SPECIES - 2]);               // H2
    xi[1] = _DF(0.15) * (xrest - xi[NUM_SPECIES - 2]);               // O2
    xi[NUM_SPECIES - 3] = _DF(0.56) * (xrest - xi[NUM_SPECIES - 2]); // Xe
#endif
    get_yi(xi, thermal.Wi);

    if (bl.RSources)
    {
        real_t xre = _DF(1.0e-15), ratios = xre * real_t(NUM_COP - 3) * _DF(0.25);
        for (size_t n1 = 0; n1 < NUM_COP; n1++)
            xi[n1] -= ratios;
        for (size_t nn = 2; nn < NUM_COP - 2; nn++)
            xi[nn] = xre;
    }

#endif // end COP
}

/**
 * @brief  Initialize conservative quantity;
 * @return void
 */
extern void InitialUFKernel(int i, int j, int k, Block bl, MaterialProperty material, Thermal thermal, real_t *U, real_t *U1, real_t *LU,
                            real_t *FluxF, real_t *FluxG, real_t *FluxH, real_t *FluxFw, real_t *FluxGw, real_t *FluxHw,
                            real_t *u, real_t *v, real_t *w, real_t *rho, real_t *p, real_t *_y, real_t *T, real_t *H, real_t *c)
{
    MARCO_DOMAIN_GHOST();
    if (i >= Xmax)
        return;
    if (j >= Ymax)
        return;
    if (k >= Zmax)
        return;

    int id = Xmax * Ymax * k + Xmax * j + i;
    real_t dx = bl.dx, dy = bl.dy, dz = bl.dz;
    real_t x = bl.DimX ? (i - Bwidth_X + bl.myMpiPos_x * (Xmax - Bwidth_X - Bwidth_X)) * dx + _DF(0.5) * dx + bl.Domain_xmin : _DF(0.0);
    real_t y = bl.DimY ? (j - Bwidth_Y + bl.myMpiPos_y * (Ymax - Bwidth_Y - Bwidth_Y)) * dy + _DF(0.5) * dy + bl.Domain_ymin : _DF(0.0);
    real_t z = bl.DimZ ? (k - Bwidth_Z + bl.myMpiPos_z * (Zmax - Bwidth_Z - Bwidth_Z)) * dz + _DF(0.5) * dz + bl.Domain_zmin : _DF(0.0);

    // Get R of mixture
    real_t *yi = &(_y[NUM_SPECIES * id]), R = get_CopR(thermal._Wi, yi);

    rho[id] = p[id] / R / T[id];

    real_t Gamma_m = get_CopGamma(thermal, yi, T[id]);
    c[id] = sycl::sqrt(p[id] / rho[id] * Gamma_m);

    // U[4] of mixture differ from pure gas
    real_t h = get_Coph(thermal, yi, T[id]);
    U[Emax * id + 4] = rho[id] * (h + _DF(0.5) * (u[id] * u[id] + v[id] * v[id] + w[id] * w[id])) - p[id];
    H[id] = (U[Emax * id + 4] + p[id]) / rho[id];
    // initial U[0]-U[3]
    U[Emax * id + 0] = rho[id];
    U[Emax * id + 1] = rho[id] * u[id];
    U[Emax * id + 2] = rho[id] * v[id];
    U[Emax * id + 3] = rho[id] * w[id];

    // initial flux terms F, G, H
    FluxF[Emax * id + 0] = U[Emax * id + 1];
    FluxF[Emax * id + 1] = U[Emax * id + 1] * u[id] + p[id];
    FluxF[Emax * id + 2] = U[Emax * id + 1] * v[id];
    FluxF[Emax * id + 3] = U[Emax * id + 1] * w[id];
    FluxF[Emax * id + 4] = (U[Emax * id + 4] + p[id]) * u[id];

    FluxG[Emax * id + 0] = U[Emax * id + 2];
    FluxG[Emax * id + 1] = U[Emax * id + 2] * u[id];
    FluxG[Emax * id + 2] = U[Emax * id + 2] * v[id] + p[id];
    FluxG[Emax * id + 3] = U[Emax * id + 2] * w[id];
    FluxG[Emax * id + 4] = (U[Emax * id + 4] + p[id]) * v[id];

    FluxH[Emax * id + 0] = U[Emax * id + 3];
    FluxH[Emax * id + 1] = U[Emax * id + 3] * u[id];
    FluxH[Emax * id + 2] = U[Emax * id + 3] * v[id];
    FluxH[Emax * id + 3] = U[Emax * id + 3] * w[id] + p[id];
    FluxH[Emax * id + 4] = (U[Emax * id + 4] + p[id]) * w[id];

#ifdef COP
    for (size_t ii = 5; ii < Emax; ii++)
    { // equations of species
        U[Emax * id + ii] = rho[id] * yi[ii - 5];
        FluxF[Emax * id + ii] = rho[id] * u[id] * yi[ii - 5];
        FluxG[Emax * id + ii] = rho[id] * v[id] * yi[ii - 5];
        FluxH[Emax * id + ii] = rho[id] * w[id] * yi[ii - 5];
    }
#endif // end COP

    // give intial value for the interval matrixes
    for (int n = 0; n < Emax; n++)
    {
        LU[Emax * id + n] = _DF(0.0);         // incremental of one time step
        U1[Emax * id + n] = U[Emax * id + n]; // intermediate conwervatives
        FluxFw[Emax * id + n] = _DF(0.0);     // numerical flux F
        FluxGw[Emax * id + n] = _DF(0.0);     // numerical flux G
        FluxHw[Emax * id + n] = _DF(0.0);     // numerical flux H
    }
}
